set.seed(5168)
## responsibility attribution ##

	library(foreign)
	library(mnormt)
	library(ordinal)

  data <- read.dta('Denmark (2014).dta')

  data <- subset(data, data$Party != 1) #Remove Venstre
  data <- subset(data, data$Party != 2) #Remove Socialdemokraterne
  data <- subset(data, data$Party != 3) #Remove Dansk Folkeparti
  data <- subset(data, data$Party != 4) #Remove Radikale
  data <- subset(data, data$Party != 7) #Remove Liberal Alliance
  data <- subset(data, data$Party != 8) #Remove Det Konservative Folkeparti
  data <- subset(data, data$Party != 9) #Remove Kristendemokraterne

	model <- clmm(as.factor(Resp_attribution) ~ Prime_minister + Cabinet + Opposition + DK_Perc_seats + Prime_minister * DK_Perc_seats + Cabinet * DK_Perc_seats + Opposition*DK_Perc_seats + (1 | pid), data = data, HESS = TRUE, link = 'probit')
	summary(model)

	pars <- as.matrix(model$coefficients)
	vce <- vcov(model)
	vce <- vce[-c(12), ]       
	vce <- vce[, -c(12)]       
	
	simCoef <- rmnorm(1000, as.vector(pars),as.matrix(vce))

	summary(model)
	summary(simCoef)

	pmHi <- c()
	pmLo <- c()
	pmMed <- c()

  partnerHi <- c()
  partnerLo <- c()
  partnerMed <- c()

	supHi <- c()
	supLo <- c()
	supMed <- c()

	oppHi <- c()
	oppLo <- c()
	oppMed <- c()

  otherHi <- c()
  otherLo <- c()
  otherMed <- c()


	for(i in 0:50){
	## Prime minister 

		agg <- simCoef[, 5] * 1 + 
				simCoef[, 6] * 0 + 
				simCoef[, 7] * 0 + 
				simCoef[, 8] * i + 
				simCoef[, 9] * i + 
				simCoef[, 10] * 0 + 
				simCoef[, 11] * 0 
		
		pred1 <- pnorm(simCoef[, 1], agg, sd = 1)
		pred2 <- pnorm(simCoef[, 2], agg, sd = 1) - pred1
		pred3 <- pnorm(simCoef[, 3], agg, sd = 1) - pred2 - pred1
		pred4 <- pnorm(simCoef[, 4], agg, sd = 1) - pred3 - pred2 - pred1
		pred5 <- 1 - pred1 - pred2 - pred3 - pred4
		
		pred <- pred5 * 5 + 
			pred4 * 4 + 
			pred3 * 3 + 
			pred2 * 2 + 
			pred1

		pmHi[i + 1] <- quantile(pred, 0.975)
		pmLo[i + 1] <- quantile(pred, 0.025)
		pmMed[i + 1] <- quantile(pred, 0.5)


	## partner
	agg <- simCoef[, 5] * 0 + 
	  simCoef[, 6] * 1 + 
	  simCoef[, 7] * 0 + 
	  simCoef[, 8] * i + 
	  simCoef[, 9] * 0 + 
	  simCoef[, 10] * i + 
	  simCoef[, 11] * 0  
		
		pred1 <- pnorm(simCoef[, 1], agg, sd = 1)
		pred2 <- pnorm(simCoef[, 2], agg, sd = 1) - pred1
		pred3 <- pnorm(simCoef[, 3], agg, sd = 1) - pred2 - pred1
		pred4 <- pnorm(simCoef[, 4], agg, sd = 1) - pred3 - pred2 - pred1
		pred5 <- 1 - pred1 - pred2 - pred3 - pred4
		
		pred <- pred5 * 5 + 
			pred4 * 4 + 
			pred3 * 3 + 
			pred2 * 2 + 
			pred1

		partnerHi[i + 1] <- quantile(pred, 0.975)
		partnerLo[i + 1] <- quantile(pred, 0.025)
		partnerMed[i + 1] <- quantile(pred, 0.5)

	## opposition
	agg <- simCoef[, 5] * 0 + 
	  simCoef[, 6] * 0 + 
	  simCoef[, 7] * 1 + 
	  simCoef[, 8] * i + 
	  simCoef[, 9] * 0 + 
	  simCoef[, 10] * 0 + 
	  simCoef[, 11] * i 
		
		pred1 <- pnorm(simCoef[, 1], agg, sd = 1)
		pred2 <- pnorm(simCoef[, 2], agg, sd = 1) - pred1
		pred3 <- pnorm(simCoef[, 3], agg, sd = 1) - pred2 - pred1
		pred4 <- pnorm(simCoef[, 4], agg, sd = 1) - pred3 - pred2 - pred1
		pred5 <- 1 - pred1 - pred2 - pred3 - pred4
		
		pred <- pred5 * 5 + 
			pred4 * 4 + 
			pred3 * 3 + 
			pred2 * 2 + 
			pred1

		oppHi[i + 1] <- quantile(pred, 0.975)
	  oppLo[i + 1] <- quantile(pred, 0.025)
	  oppMed[i + 1] <- quantile(pred, 0.5)

    ## No seats
	agg <- simCoef[, 5] * 0 + 
	  simCoef[, 6] * 0 + 
	  simCoef[, 7] * 0 + 
	  simCoef[, 8] * i + 
	  simCoef[, 9] * 0 + 
	  simCoef[, 10] * 0 + 
	  simCoef[, 11] * 0 
    
    pred1 <- pnorm(simCoef[, 1], agg, sd = 1)
    pred2 <- pnorm(simCoef[, 2], agg, sd = 1) - pred1
    pred3 <- pnorm(simCoef[, 3], agg, sd = 1) - pred2 - pred1
    pred4 <- pnorm(simCoef[, 4], agg, sd = 1) - pred3 - pred2 - pred1
    pred5 <- 1 - pred1 - pred2 - pred3 - pred4
    
    pred <- pred5 * 5 + 
      pred4 * 4 + 
      pred3 * 3 + 
      pred2 * 2 + 
      pred1
    
    otherHi[i + 1] <- quantile(pred, 0.975)
    otherLo[i + 1] <- quantile(pred, 0.025)
    otherMed[i + 1] <- quantile(pred, 0.5)
    
    }


## now paint a pretty picture.... ##

	seats <- seq(0, 50, 1)        
  	
	#Now create the figure for Denmark
	par(mar=c(3.6,4,1,2))      
	
	plot(1, 1, type = 'n',
		xlim = c(0, 50),
		ylim = c(1, 5),
		xlab = '',
		ylab = '',
		main = '',
		axes = FALSE)
	axis(1,at=seq(0,50,25),labels=T)         
	axis(2)

  lines(seats, partnerMed, col = rgb(0, 0, 1, 0.2), lwd = 3)

	seatHi <- quantile(data$DK_Perc_seats[data$Cabinet == 1], 0.9, na.rm = TRUE)
	seatLo <- quantile(data$DK_Perc_seats[data$Cabinet == 1], 0.1, na.rm = TRUE)
	lines(seats[seats >= seatLo & seats <= seatHi], partnerMed[seats >= seatLo & seats <= seatHi], col = rgb(0, 0, 1, 0.3), lwd = 15)

	lines(seats, oppMed, col = rgb(0, 0, 0, 0.2), lwd = 3)

	seatHi <- quantile(data$DK_Perc_seats[data$Opposition == 1], 0.9, na.rm = TRUE)
	seatLo <- quantile(data$DK_Perc_seats[data$Opposition == 1], 0.1, na.rm = TRUE)
	lines(seats[seats >= seatLo & seats <= seatHi], oppMed[seats >= seatLo & seats <= seatHi], col = rgb(0, 0, 0, 0.3), lwd = 15)

	mtext("Responsibility attribution", side=2, line=2.5, cex=1.5)     
	mtext("Perceived seat %", side=1, line=2.2, cex=1.5)               


############
###LEGEND###
############

legend('bottomright', legend = c('Partner', 'Opposition'), fill = c(rgb(0, 0, 1, 0.3), rgb(0, 0, 0, 0.3)), bty = 'n', border = NA, , cex=1.25)

###################
###EXPORT FIGURE###
###################

dev.copy(png,'Figure 9b DK 2014 - Effect of support party categorization.png', height = 490, width = 633)	
dev.off()  
